Intro a series de tiempo

Author

Verónica Andreo

Published

August 21, 2023

Procesamiento de series de tiempo en GRASS GIS

Contenidos

  • Nociones básicas
  • TGRASS framework
  • Creación de series de tiempo
  • Álgebra temporal y variables temporales
  • Distintos tipos de agregación
  • Estadística zonal e islas de calor urbanas
  • Conexión con R

GRASS GIS es el primer SIG de código abierto que incorporó capacidades para gestionar, analizar, procesar y visualizar datos espacio-temporales, así como las relaciones temporales entre series de tiempo.

TGRASS: GRASS Temporal

  • Completamente basado en metadatos, por lo que no hay duplicación de datos
  • Sigue una aproximación Snapshot, i.e., añade marcas de tiempo o timestamps a los mapas
  • Una colección de mapas de la misma variable con timestamps se llama space-time dataset o STDS
  • Los mapas en una STDS pueden tener diferentes extensiones espaciales y temporales
  • TGRASS utiliza una base de datos SQLite para almacenar la extensión temporal y espacial de las STDS, así como las relaciones topológicas entre los mapas y entre las STDS en cada mapset.

Space-time datasets

  • Space time raster datasets (STRDS)
  • Space time 3D raster datasets (STR3DS)
  • Space time vector datasets (STVDS)

Otras nociones básicas en TGRASS

  • El tiempo puede definirse como intervalos (inicio y fin) o como instancias (sólo inicio)
  • El tiempo puede ser absoluto (por ejemplo, 2017-04-06 22:39:49) o relativo (por ejemplo, 4 años, 90 días)
  • Granularidad es el mayor divisor común de todas las extensiones temporales (y posibles gaps) de los mapas de un STDS
  • Topología se refiere a las relaciones temporales entre los intervalos de tiempo en una STDS

  • Muestreo temporal se utiliza para determinar el estado de un proceso durante un segundo proceso.

Módulos temporales

  • t.*: Módulos generales para manejar STDS de todos los tipos
  • t.rast.*: Módulos que tratan con STRDS
  • t.rast3d.*: Módulos que tratan con STR3DS
  • t.vect.*: Módulos que tratan con STVDS

TGRASS: marco general y flujo de trabajo

Manos a la obra con series de tiempo raster en GRASS GIS

Datos para la sesión

  • Producto MODIS: MOD11B3 Collection 6
  • Tile: h12v12
  • Composiciones mensuales
  • Resolución espacial: 5600m
  • Mapset modis_lst

Código para la sesión

Iniciar GRASS GIS directamente en el mapset modis_lst

# paths
grassdata='/home/veroandreo/grassdata/'
location='posgar2007_4_cba'
mapset='modis_lst'
import os
import subprocess
import sys

# Ask GRASS GIS where its Python packages are to be able to start it from the notebook
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

# Importar los paquetes python de GRASS
import grass.script as gs
import grass.jupyter as gj

# Iniciar GRASS
session = gj.init(grassdata, location, mapset)

Establecer región computacional y máscara

Listar los mapas raster y obtener información de uno de ellos

# Get list of raster maps in the 'modis_lst' mapset
g.list type=raster mapset=.

# Get info from one of the raster maps
r.info map=MOD11B3.A2015060.h12v12.single_LST_Day_6km

Establecer la región computacional

# Set region to Cba boundaries with LST maps' resolution
g.region -p vector=provincia_cba \
  align=MOD11B3.A2015060.h12v12.single_LST_Day_6km

Aplicar máscara

# Set a MASK to Cba boundary
r.mask vector=provincia_cba

Crear un conjunto de datos espacio-temporales (STDS)

t.create
- Crea una tabla SQLite en la base de datos temporal - Permite manejar grandes cantidades de mapas usando el STDS como entrada - Necesitamos especificar: - tipo de mapas (raster, raster3d o vector) - tipo de tiempo (absoluto o relativo)

Crear la STRDS

# Create the STRDS
t.create type=strds temporaltype=absolute output=LST_Day_monthly \
  title="Monthly LST Day 5.6 km" \
  description="Monthly LST Day 5.6 km MOD11B3.006 Cordoba, 2015-2019"

Chequear si la STRDS fue creada

# Check if the STRDS is created
t.list type=strds

Obtener información sobre la STRDS

# Get info about the STRDS
t.info input=LST_Day_monthly

Registrar mapas en una STDS (asignar timestamps)

t.register
- Asigna o agrega timestamps a los mapas - Necesitamos: - el STDS vacío como entrada, i.e., la tabla SQLite contenedora, - la lista de mapas que se registrarán, - la fecha de inicio, - la opción de incremento junto con -i para la creación de intervalos

Añadir timestamps a los mapas, i.e., registrar mapas

# Add time stamps to maps (i.e., register maps)
t.register -i input=LST_Day_monthly \
 maps=`g.list type=raster pattern="MOD11B3*LST_Day*" separator=comma` \
 start="2015-01-01" increment="1 months"

Chequear la información sobre la STRDS nuevamente

# Check info again
t.info input=LST_Day_monthly

Obtener la lista de mapas en la STRDS

# Check the list of maps in the STRDS
t.rast.list input=LST_Day_monthly

Chequear los valores mínimos y máximos de cada mapa

# Check min and max per map
t.rast.list input=LST_Day_monthly columns=name,min,max

Para más opciones, ver el manual de t.register y la wiki sobre opciones para registrar mapas en STDS.

Representación gráfica de STDS

Crear una representación gráfica de la serie de tiempo

# graphical representation of our STRDS
g.gui.timeline inputs=LST_Day_monthly

Ver el manual de g.gui.timeline

Operaciones con álgebra temporal: t.rast.algebra

  • Realiza una amplia gama de operaciones de álgebra temporal y espacial basadas en la topología temporal de los mapas
    • Operadores temporales: unión, intersección, etc.
    • Funciones temporales: start_time(), start_doy(), etc.
    • Operadores espaciales (subconjunto de r.mapcalc)
    • Modificador de vecindario temporal: [x,y,t]
    • Otras funciones temporales como t_snap(), buff_t() o t_shift()

¡pueden combinarse en expresiones complejas!

Desde K*50 a Celsius usando la calculadora temporal

Re-escalar a grados Celsius

# Re-scale data to degrees Celsius
t.rast.algebra basename=LST_Day_monthly_celsius suffix=gran \
  expression="LST_Day_monthly_celsius = LST_Day_monthly * 0.02 - 273.15"

Ver info de la nueva serie de tiempo

# Check info
t.info LST_Day_monthly_celsius

Gráfico temporal: LST vs tiempo

Gráfico temporal de LST para la ciudad de Córdoba, Argentina

# LST time series plot for Cba city center
g.gui.tplot strds=LST_Day_monthly_celsius \
  coordinates=4323478.531282977,6541664.09350761 \
  title="Monthly LST. City center of Cordoba" \
  xlabel="Time" ylabel="LST"

Para un único punto, ver g.gui.tplot. Para un vector de puntos, ver t.rast.what.

Las coordenadas del punto pueden ser escritas directamente, copiadas desde el mapa o seleccionadas interactivamente.

Listas y selecciones

  • t.list para listar las STDS y los mapas registrados en la base de datos temporal,
  • t.rast.list para mapas en series temporales de rasters, y
  • t.vect.list para mapas en series temporales de vectores.

Variables usadas para hacer las listas y selecciones

STRDS:id, name, creator, mapset, temporal_type, creation_time, start_time, end_time, north, south, west, east, nsres, ewres, cols, rows, number_of_cells, min, max

STVDS:id, name, layer, creator, mapset, temporal_type, creation_time, start_time, end_time, north, south, west, east, points, lines, boundaries, centroids, faces, kernels, primitives, nodes, areas, islands, holes, volumes

Ejemplos de listas y selecciones

Mapas cuyo valor mínimo es menor o igual a 10

# Maps with minimum value lower than or equal to 10
t.rast.list input=LST_Day_monthly_celsius order=min \
 columns=name,start_time,min where="min <= '10.0'"

Mapas cuyo valor máximo es mayor a 30

# Maps with maximum value higher than 30
t.rast.list input=LST_Day_monthly_celsius order=max \
 columns=name,start_time,max where="max > '30.0'"

Mapas contenidos entre dos fechas

# Maps between two given dates
t.rast.list input=LST_Day_monthly_celsius columns=name,start_time \
 where="start_time >= '2015-05' and start_time <= '2015-08-01 00:00:00'"

Todos los mapas correspondientes al mes de Enero

# Maps from January
t.rast.list input=LST_Day_monthly_celsius columns=name,start_time \
 where="strftime('%m', start_time)='01'"

Estadística descriptiva de STRDS

Imprimir estadísticas descriptivas univariadas para cada mapa dentro de la STRDS

# Print univariate stats for maps within STRDS
t.rast.univar input=LST_Day_monthly_celsius

Obtener estadísticas extendidas con la opción -e

# Get extended statistics
t.rast.univar -e input=LST_Day_monthly_celsius

Escribir la salida a un archivo de texto

# Write the univariate stats output to a csv file
t.rast.univar input=LST_Day_monthly_celsius separator=comma \
  output=stats_LST_Day_monthly_celsius.csv

Agregación temporal 1: Serie completa

t.rast.series

  • Agrega STRDS completas o partes de ellas usando la opción where.
  • Diferentes métodos disponibles: promedio, mínimo, máximo, mediana, moda, etc.

LST máxima y mínima del período 2015-2019

Obtener el mapa de la máxima LST del período

# Get maximum LST in the STRDS
t.rast.series input=LST_Day_monthly_celsius \
  output=LST_Day_max method=maximum

Obtener el mapa de la mínima LST del período

# Get minimum LST in the STRDS
t.rast.series input=LST_Day_monthly_celsius \
  output=LST_Day_min method=minimum

Cambiar la paleta de colores a celsius

# Change color pallete to celsius
r.colors map=LST_Day_min,LST_Day_max color=celsius

Comparar mapas con la herramienta Mapswipe

Operaciones usando variables temporales

t.rast.mapcalc

  • Ejecuta expresiones espacio-temporales tipo r.mapcalc
  • Permite operadores espaciales y temporales, así como variables internas en la expresión
  • Las variables temporales incluyen: start_time(), end_time(), start_month(), start_doy(), etc.

Cuál es el mes de máxima LST?

Obtener el mes en que ocurre el máximo de LST en cada pixel

# Get month of maximum LST
t.rast.mapcalc -n inputs=LST_Day_monthly_celsius \
  output=month_max_lst \
  expression="if(LST_Day_monthly_celsius == LST_Day_max, start_month(), null())" \
  basename=month_max_lst

Obtener información del mapa resultante

# Get basic info
t.info month_max_lst

Obtener el primer mes en que aparece el máximo de LST

# Get the earliest month in which the maximum appeared (method minimum)
t.rast.series input=month_max_lst \
  method=minimum \
  output=max_lst_date

Remover la STRDS intermedia y los mapas que contiene: month_max_lst

# Remove month_max_lst strds 
# we were only interested in the resulting aggregated map
t.remove -rf inputs=month_max_lst

Mostrar el mapa resultante desde la terminal

Abrir un monitor wx

# Open a monitor
d.mon wx0

Mostrar el mapa raster

# Display the raster map
d.rast map=max_lst_date

Mostrar sólo los bordes del mapa vectorial de NC

# Display boundary vector map
d.vect map=provincia_cba type=boundary color=#4D4D4D width=2

Agregar leyenda

# Add raster legend
d.legend -t raster=max_lst_date title="Month" \
  labelnum=6 title_fontsize=20 font=sans fontsize=16

Agregar barra de escala

# Add scale bar
d.barscale length=100 units=kilometers segment=4 fontsize=14

Agregar Norte

# Add North arrow
d.northarrow style=1b text_color=black

Agregar título

# Add text
d.text text="Month of maximum LST" \
  color=black align=cc font=sans size=12

Podríamos haber hecho lo mismo pero anualmente para conocer en qué mes ocurre el máximo en cada año y así evaluar la ocurrencia de tendencias. Cómo lo harían?

Agregación temporal 2: granularidad

t.rast.aggregate

  • Agrega mapas raster dentro de STRDS con diferentes granularidades
  • La opción where permite establecer fechas específicas para la agregación
  • Diferentes métodos disponibles: promedio, mínimo, máximo, mediana, moda, etc.

De LST mensual a estacional

LST media estacional

# 3-month mean LST
t.rast.aggregate input=LST_Day_monthly_celsius \
  output=LST_Day_mean_3month \
  basename=LST_Day_mean_3month suffix=gran \
  method=average granularity="3 months"

Chequear info

# Check info
t.info input=LST_Day_mean_3month

Chequear lista de mapas

# Check map list
t.rast.list input=LST_Day_mean_3month

Tarea

Comparar las líneas de tiempo mensual y estacional con g.gui.timeline

g.gui.timeline inputs=LST_Day_monthly_celsius,LST_Day_mean_3month

Graficar LST estacional con monitores wx

Establecer la paleta de colores celsius para la STRDS estacional

# Set STRDS color table to celsius degrees
t.rast.colors input=LST_Day_mean_3month color=celsius

Iniciar un monitor Cairo

# Start a new graphics monitor
d.mon cairo out=frames.png width=1400 height=500 resolution=4 --o

Crear el primer frame

# create a first frame
d.frame -c frame=first at=0,100,0,25
d.rast map=LST_Day_mean_3month_2015_01
d.vect map=provincia_cba type=boundary color=#4D4D4D width=2
d.text text='Ene-Mar 2015' color=black font=sans size=6 bgcolor=white

Crear el segundo frame

# create a second frame
d.frame -c frame=second at=0,100,25,50
d.rast map=LST_Day_mean_3month_2015_04
d.vect map=provincia_cba type=boundary color=#4D4D4D width=2
d.text text='Abr-Jun 2015' color=black font=sans size=6 bgcolor=white

Crear el tercer frame

# create a third frame
d.frame -c frame=third at=0,100,50,75
d.rast map=LST_Day_mean_3month_2015_07
d.vect map=provincia_cba type=boundary color=#4D4D4D width=2
d.text text='Jul-Sep 2015' color=black font=sans size=6 bgcolor=white

Crear el cuarto frame

# create a fourth frame
d.frame -c frame=fourth at=0,100,75,100
d.rast map=LST_Day_mean_3month_2015_10
d.vect map=provincia_cba type=boundary color=#4D4D4D width=2
d.text text='Oct-Dic 2015' color=black font=sans size=6 bgcolor=white

Liberar el monitor

# release monitor
d.mon -r

LST estacional en 2015

Tarea

Ahora que ya conocen t.rast.aggregate, extraigan el mes de máximo LST por año y luego vean si hay alguna tendencia positiva o negativa, es decir, si los valores máximos de LST se observan más tarde o más temprano con el tiempo (años)

Una solución podría ser…

t.rast.aggregate \
  input=LST_Day_monthly_celsius \
  output=month_max_LST_per_year \
  basename=month_max_LST suffix=gran \
  method=max_raster \
  granularity="1 year" 

t.rast.series \
  input=month_max_LST_per_year \
  output=slope_month_max_LST \
  method=slope

Animaciones

Animación de la serie estacional de LST

# Animation of seasonal LST
g.gui.animation strds=LST_Day_mean_3month

Ver el manual de g.gui.animation para más opciones y ajustes.

Agregación vs Climatología

Agregación por granularidad

Agregación tipo climatología

Climatologías mensuales

LST promedio de Enero

# January average LST
t.rast.series input=LST_Day_monthly_celsius \
  method=average \
  where="strftime('%m', start_time)='01'" \
  output=LST_average_jan

Climatología para todos los meses

# for all months - *nix
for MONTH in `seq -w 1 12` ; do 
 t.rast.series input=LST_Day_monthly_celsius method=average \
  where="strftime('%m', start_time)='${MONTH}'" \
  output=LST_average_${MONTH}
done

Tarea

  • Comparar las medias mensuales con las climatologías mensuales
  • Las climatologías que creamos forman una STRDS?

Anomalías anuales

\[AnomaliaStd_i = \frac{Media_i - Media}{SD}\]

Se necesitan:

  • promedio y desviación estándar general de la serie
  • promedios anuales

Obtener el promedio general de la serie

# Get general average
t.rast.series input=LST_Day_monthly_celsius \
 method=average output=LST_average

Obtener el desvío estándar general de la serie

# Get general SD
t.rast.series input=LST_Day_monthly_celsius \
 method=stddev output=LST_sd

Obtener los promedios anuales

# Get annual averages
t.rast.aggregate input=LST_Day_monthly_celsius \
 method=average granularity="1 years" \
 output=LST_yearly_average basename=LST_yearly_average

Estimar las anomalías anuales

# Estimate annual anomalies
t.rast.algebra basename=LST_year_anomaly \
 expression="LST_year_anomaly = (LST_yearly_average - map(LST_average)) / map(LST_sd)"

Establecer la paleta de colores differences

# Set difference color table
t.rast.colors input=LST_year_anomaly color=difference

Animación

# Animation of annual anomalies
g.gui.animation strds=LST_year_anomaly

Isla de calor superficial urbana (Surface Urban Heat Island - SUHI)

  • La temperatura del aire de una zona urbana es más alta que la de las zonas cercanas
  • La UHI tiene efectos negativos en la calidad del agua y el aire, la biodiversidad, la salud humana y el clima.
  • La SUHI también está muy relacionada con la salud, ya que influye en la UHI

SUHI y área rural en Buenos Aires. Fuente: Wu et al, 2019.

Estadística zonal en series de tiempo de datos raster

v.strds.stats - Permite obtener datos de series de tiempo agregados espacialmente para polígonos de un mapa vectorial

SUHI estival para Córdoba y alrededores

Instalar la extensión v.strds.stats

# Install v.strds.stats add-on
g.extension extension=v.strds.stats

Listar mapas

# List maps in seasonal time series
t.rast.list input=LST_Day_mean_3month

Extraer LST promedio de verano para el Gran Córdoba

# Extract summer average LST for Cba urban area
v.strds.stats input=area_edificada_cba \
  strds=LST_Day_mean_3month \
  where="fna == 'Gran Córdoba'" \
  t_where="strftime('%m', start_time)='01'" \
  output=cba_summer_lst \
  method=average

Crear buffer externo - 30 km

# Create outside buffer - 30 km
v.buffer input=cba_summer_lst \
  distance=30000 \
  output=cba_summer_lst_buf30

Crear buffer interno - 15 km

# Create inside buffer - 15 km
v.buffer input=cba_summer_lst \
  distance=15000 \
  output=cba_summer_lst_buf15

Remover el área del buffer 15 km del buffer de 30 km

# Remove 15km buffer area from the 30km buffer area
v.overlay ainput=cba_summer_lst_buf15 \
  binput=cba_summer_lst_buf30 \
  operator=xor \
  output=cba_surr

Límites del Gran Córdoba y el área rural circundante

Extraer estadísticas para los alrededores del Gran Córdoba

# Extract zonal stats for Cba surroundings
v.strds.stats input=cba_surr \
  strds=LST_Day_mean_3month \
  t_where="strftime('%m', start_time)='01'" \
  method=average \
  output=cba_surr_summer_lst

Chequear la LST estival promedio para el Gran Córdoba y alrededores

# Take a look at mean summer LST in Cba and surroundings
v.db.select cba_summer_lst
v.db.select cba_surr_summer_lst

Recursos (muy) útiles

Referencias

  • Gebbert, S., Pebesma, E. (2014) A temporal GIS for field based environmental modeling. Environmental Modelling & Software, 53, 1-12. DOI
  • Gebbert, S., Pebesma, E. (2017) The GRASS GIS temporal framework. IJGIS 31, 1273-1292. DOI
  • Gebbert, S., Leppelt, T., Pebesma, E. (2019) A Topology Based Spatio-Temporal Map Algebra for Big Data Analysis. Data, 4, 86. DOI